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ABSTRACT 

In  this  paper  we  show  how  to  modify  the  random  choice 
method  of  Glimm  by  replacing  the  exact  solution  of  the  Riemann 
problem  with  an  appropriate  finite  difference  approximation. 
Our  modification  resolves  discontinuities  as  well  as  Glimm' s 
scheme,  but  is  computationally  more  efficient  and  is  easier 
to  extend  to  more  general  situations. 


Introduction 

Glimm's  random  choice  method  is  a  semi-analytical 
technique  for  calculating  discontinuous  solutions  of  hyperbolic 
systems  of  conservation  laws.   Approximate  solutions  are 
represented  as  piecewise  constant  at  any  time;  they  are 
advanced  in  time  by  solving  exactly  the  Riemann  initial 
value  problem   formed  by  the  constant  states  between  two 
neighboring  cells.   The  value  of  the  approximation   in  each 
cell  at  the  new  time   is  taken  to  be  the  value  of  the  exact 
solution  at  a  randomly  chosen  point  in  the  cell.   Glimm  has 
proved  the  convergence  of  the  method  in  one  space  dimension 
and  for  initial  data  which  differ  little  from  a  constant  state. 

The  original  version  of  Glimm's  method  had  an  overall 
accuracy  of  order  half;  this  was  later  improved  by  Glimm  to 
accuracy  of  first  order.   The  main  advantages  of  the  method 
for  numerical  calculations  are  the  sharp  resolution   of  dis- 
continuities and  absence  of  over   and  undershoots.   Chorin 
has  pointed  out  that  these  features  were  decisive  for  calcula- 
tions of  combustion  problems,  and  has  made  significant  use 
of  the  method  for  these  problems,  [1],  [2]. 

A  drawback  of  the  method  is  the  need  to  solve  exactly 
many  Riemann  problems.   This  is  a  time  consuming  procedure 
and  furthermore  stands  in  the  way  of  extending  the  method  to 
situations  where  exact  solutions  are  not  available.  The  goal 
of  this  paper  is  to  show  how  to  replace  the  exact  solution  of 
the  Riemann  problem  by  an  approximate  one  while  retaining  the 
advantages  of  Glimm's  method. 


In  the  course  of  analyzing  our  method  we  estimate 
the  probability  of  gross  error   by  the  random  choice  method; 
this  estimate  is  a  great  improvement  of  the  (alarmingly) 
high  previous  estimates. 

The   organization  of  the  paper  is  as  follows: 

1.  Review  of  theory  of  weak  solutions. 

2.  Glimm's  scheme  and  its  modifications. 

3.  Convergence  of  the  random  choice  method. 

4.  A  finite  difference  approximation  of  the  Riemann 
problem. 

5.  Scalar  conservation  laws. 

In  a  sequel  we  shall  apply  these  ideas  to  the  equations 
of  compressible  flow. 


1 .     Introduction.   Review  of  Theory  of  Weak  Solutions. 

In  this  paper  we  consider  numerical  solutions  of  the 
initial   value  problem  for  hyperbolic  systems  of  conservation 
laws . 


(1.1a)   u   +  f(u)^  =  0  ,    u(x,0)  =  ())(x)  ,    -co  <  X  <  «>  . 


Here  u(x,t)  is  an  m-column  vector  of  unknowns  and  f(u),  the 
flux,  is  a  vector  valued  function  of  m  components.  We  can 
write  (1.1)   in  matrix  form. 

(1.1b)  u^  +  Au   =  0  , 

t      X      ' 

where 
(1.2)  A(u)  =  f^   . 

(1.1)  is  called  hyperbolic  if  all  eigenvalues  of  the  matrix  A 
are  real. 

It  is  well  known  that  solutions  to  (1.1)   may  develop 
shocks  and  contact  discontinuities  even  when  the  initial  data 
is  smooth.   To  allow  for  discontinuous  solutions   one  admits 
weak  solutions   which  satisfy  the  system  (1.1)   in  the  sense 
of  distribution  theory,  i.e. 


(1.3)        [w  u  +  w^f(u)J  dx  dt  +     w(x,0)  u(x,0)  dx  =  0 

0   -co  _oo 

for  all  C   that  functions  w(x,t)   test  vanishes  for  | x | +t  large, 
Clearly,  a  piecewise  smooth  u  is  a  wear  -solution  of  (1.1)  if 
and  only  if 


(i)    u  satisfies  u   +  f(u)   =  0   pointwise  in  each  smooth  region, 
(ii)   across  each  curve  of  discontinuity   the  Rankine-Hugoniot 
relation  (abbreviated  RH) 

(1.4)  f(u^)  -  f(Uj^)  =  S(Uj^  -  u^) 

holds,  where  S  is  the  speed  of  propagation  of  the  discontinuity, 

and  u   and  u   are  the  states  on  the  left  and  the  right  of  the 
L        R 

discontinuity,  respectively. 

Weak  solutions  of  (1.1)   are  not  uniquely  determined  by 
their  initial  data;  an  additional  principal  is  needed  to  select 
the  physically  relevant  solution.   One  such  selection  principle 
is  to  admit  only  those  weak  solutions  u  which  are  limits  as 
e  ^  0  of  solutuons   u(e)  of  the  viscous  equations 


(1.5)  u^  +  f(u)^  =   £U^^  ,  £  >  0. 


Such  limit  solutions  can  be  characterized  with  the 
aid  of  an  entropy    function    U(u)  for  (1.1)  ,  defined  as  follows: 
(i)   U  satisfies 

(1.6)  U  f   =  F 

^    '  u  u     u 

where  F  is  some  other  function  called  entropy   flux. 

(ii)   U  is  a  convex  function  of  u. 

It  follows  from  (1.6),  upon  multiplication  of  (1.1b)  by 

U   ,  that  every  smooth  solution  of  (1.1a)  also  satisfies 
u  -^ 

(1.7)  U.  +  F   =0. 

-1  t      X 


It  was  shown  in  [lol  that  if  u  is  a  limit  as  e  ->-  0 
boundediy  and  a.e.,  of  a  sequence  of  solutions  u(e)  of  (1.5), 
then  the  limit  satisfies,  in  the  weak  sense    the  following 
inequality: 


{1.8a)  U(u)^  +  F(u)  <_   0 


Weak  sense  means  that  for  all  nonnegative    smooth  test 
functions  w(x,t)  of  compact  support 

OO     OO  CO 

(1.8b)   -       [w^U  +  w^F]  dx  dt  -     w(x,0)  U(u(x,0))  dx  £ 

0   -OO  _oo 

If  u  is  piecewise  smooth  with  discontinuities,  then  (1.7) 
holds  pointwise  in  the  smooth  regions,  while  across  a 
discontinuity 


1.8c)    tF(u^)  -  F(Uj^)]  -  S[U(Uj^)  -  U(u^)]  <  0 


Relations  (1.8)  are  called  entropy    conditions . 

(1.5)  is  a  system  of  m  equations  for  two  functions 
U  and  F.  For  m  <^  2   there  are  plenty  of  solutions.  For  m  >  2 
the  existence  of  an  entropy  function  is  a  special  property 
of  the  system  (1.1).    Godunov  [7],  Friedrichs-Lax  [  5]  and 
Mock  [13]   have  observed  that  these  special  cases  can  be 
characterized  by  the  property  that  they  can  be  put  in 
symmetric  hyperbolic  form. 


A  system  of  equations 

(1.9a)  Pv^  +  Bv   =0 

t       X 

is  called  symmetric    hypevb  olio    if  P  and  B  are  symmetric 

matrices,  and  if  P  is  positive    definite . 

The  symmetrization  of  (1.1)   will  be  accomplished  by 

introducing  new  dependent  variables  v  in  place  of  u  by  setting 

u  =  u(v).   (1.1)  becomes  of  form  (1.9a)  with 

(1.9b)  P  =  u   ,     B  =  f   . 

V  V 

The  symmetry  of  the  matrices   u   ,  f    implies  that  both 
u  and  f  are  gradients  with  respect  to  v,  i.e.  there  exist 
functions    q(v)  and  p(v)  such  that 

(1.10a)  q^  -  u"^ 

(1.10b)  p   =  f"^ 

^v 

where  superscript  T  denotes  transpose.   The  positive  definite- 
ness  of  u   is  equivalent  to  the  convexity   of  q(v) . 

Note  that  the  convexity  of  q  implies  that  the  mapping 
V  ->  q  is  one-to-one,  so  that  (1.10)  can  be  inverted,  i.e. 
v  can  be  regarded  as  a  function  of  u. 

Theorems  (Godunov) .   Suppose   (1.1)  can  be  symmetrized 
by  introducing  new  variables   v,  i.e.  suppose  that  (1.10) 
holds,  where   q  is  a  convex  function  of  v.  Then  (1.1)  has 
an  entropy  function  given  by 

(1.11a)  U(u)  =  u'^v  -  q(v)  , 

with  entropy  flux, 

(1.11b)  F(u)  =  f'^v  -  p(v)  . 


Proof ;   Differentiate  (1.11a)  with  respect  to  u;  using 
(1.10a)  we  get 

T      T  T 

(1.12)  U   =  V   +  u  V   -  q  V   =  V   . 

Similarly,  from  (1.11b)  and  (1.10b)  we  get 

F   =vf   +fv,,  -pv   =vf 
u       u       u    '^v  u       u 

Relation  (1.6)  follows. 

To  prove  the  convexity  of  U,  we  show  that  U  is  the  Legendre 
transform  of  q: 

(1.13)  U(u)  =  max  [u'^v  -  q(v)]  . 

V 

For,  by  the  convexity  of  q,  the  right  side  has  a  unique  maximum; 
at  the  maximum  point  the  v  derivative  must  vanish;  this  gives 
relation  (1.10a).   This  proves  that  (1.13)  is  the  same  as  (1.11a). 

(1.13)  represents  U  as  the  maximum  of  linear  functions; 
this  proves  that  U  is  convex.   This  completes  the  proof  of 
Theorem  1.1. 

Conversely: 

Theorem  1 . 2  (Mock).   Suppose  U(u)  is  an  entropy  function 
for  (1.1)  ;  then 

(1.14)  v"^  =  U 

u 

symmetrizes  (1.1). 

Proof:   We  note  that  the  convexity  of  U  implies  that  the 
mapping  u  ->  U   is  one-to-one,  so  that  (1.14)  defines  u  as  a 
function  of  v. 


We  define  now  q  and  p  by 
(1.15a)  q(v)  =  v'^u  -  U(u) 

(1.15b)  p(v)  =  v'^f  -  F(u) 

where  F  is  entropy  flux.   Differentiating  (1.15a)  with  respect 
to  V,  and  using  (1.14)  gives 

T     T  T 

q   =u   -vu   -Uu   =u 
^v  V     u  V 

Similarly,  differentiating  (1.15b)  and  using  (1.6)  and  (1.14) 
gives 

p   =f   +vfu   -Fu   =f 
'^v  u  V     u  V 

These  formulas  show  that  (1.10a)  and  (1.10b)  hold;  therefore 

u   ,  f    are  symmetric.   To  show  that  u   is  positive  definite 
V     V        -^  v     "^ 

we  have  to  verify  that  q  is  convex.   This  can  be  done,  as  before, 
by  observing  that,  because  of  the  convexity  of  U,  it  follows 
from  (1.15a)  and  (1.14)  that  q  is  the  Legendre  transform  of  U. 
We  make  the  assumption  that  a  solution   of  the  conserva- 
tion laws  in  the  integral  sense  (1.3)  which  satisfies  the 
integral  form  (1.8b)  of  the  entropy  condition  is  uniquely 
determined  by  its  initial  values.   This  is  plausible,  on 
physical  grounds,  for  the  equations  of  compressible  flow 
but  has  not  yet  been  proven.   Uniqueness  theorems  based  on  an 
entropy  condition  for  simpler  systems  have  been  given  by 
Di  Perna,  [  4  ] . 


2.     Glimm's  Scheme  and  its  Modifications 

Glimin's  scheme  constructs  approximate  solutions  that 
are  piecewise  constant  at  each  t  =  t   on  intervals  of  length  2T : 

I.  =  ((j  -  i)r,  (j  +  i)r) 

At  the  boundary  point  of  two  adjacent  cells  there  is  a 
discontinuity.   As  in  an  earlier  method  of  Godunov,  this 
Riemann  initial  value  problem  is  solved  exactly;  the  size  t 
of  the  time  step  is  restricted   so  that  waves  issuing  from 
different  boundary  points  don't  interact: 

(2.1)  T  max  |a|  ^   T    ; 

here   max  |a|  is  the  largest   sound  speed  of  the  system  (1.1) 

at  time  t  .   Under  this  restriction  the  union  of  the  solutions 
n 

of  the  Riemann  problems  gives  an  exact  solution  at  time 

t  ,T  =  t   +  T.   Now  we  subdivide  the  real  axis  into  intervals, 
n+1     n 

and  in  each  interval  we  take  the  approximate  solution  to  be 

the  value  of  the   exact  solution  at  some  point  chosen  at  random 

in  that  interval.   In  Glimm's  original  scheme  there  was  a 

separate  random  choice  in  each  interval;  in  the  implementation 

by  Chorin  this  was  replaced  by  a  single  random  choice,  the 

same  in  all  intervals. 

Glimm's  original  version   employed  a  subdivision  at 

time  t  , ,  that  was  staggered  with  respect  to  the  subdivision 
n+1  ^^  ^ 

at  time  t  : 
n 


^n-Hl 


t 
n 


I . 


H — o — h 


D-l     D+J- 

In  tnis  case   I .   at  time   t  , ,   belongs  to  the  domain  of 
j  n+1        ^ 

deterrainacy  of   I   i  and  1- ,^    at  time  t   ;  Glimm's  scheme 
gives  the  following  approximation  v: 

(2.2)     '^^^'^^n+i^  "^  ^^^n'  ^L'  ^R^  ^°^      x  in  I   , 

where  x   is  a  point  chosen  at  random  on  I . .and   u^  ,  u_  denote 
n       "^  ]  '       L     R 

the   values  of  the  approximate  solution  v  at  time  t  on  I .  , , 


Glimm  has  snown  that  at  any  time  t   the  total  variation 

-'        n 

in  X  of  the  approximate  solutions  constructed  by  his  method 
is  uniformly  bounded,  no  matter  what  random  points  are  chosen. 
Using  this   he  has  shown  that  for  almost  all  choices  of  random 
points  the  scheme  converges  to  a  solution  of  the  conservation 
law   in  the  integral  sense  (1.3).   The  same  argument  can  be 
used  to  show  that  the  limit  solutions  satisfy  the  integral 
form  (1.8b)   of  the  entropy  condition. 

A  distinguishing  feature  of  Glimm's  scheme  is  the  absolute 
sharpness  with  which  discontinuities  are  resolved.  Consider 
e.g.  the  application   of  Glimm's  scheme  to  initial  data  that 
represent  a  single  shock: 


10 


u    for   X  <  0 
u(x,0)  =  ' 

u    for   X  >  0 

K 

The  exact  solution  of  this  Riemann  problem  at  time  x  is 

u    for   X  <  tS 
u(x,u^,Uj^)  =  ' 

[    u    for   X  >  tS 

where  the  speed  S  is  given  by  the  RH   condition  (1.4)  . 

Applying  Glimm's  recipe  (2.2)  once  has  the  effect  of 
shifting  the  initial  data  by  r  to  the  right  or  to  the  left, 
depending  on  whether  x,  is  <  xS  or  >  xS.   N  applications  give 

u    for   X  <  y 
v(x,t^)  =     L  N 

"r   f°^   ^  '  ^N 

where 

Yn  =  ^  I  sgn  (xS  -  x^)  . 

Since  each  x   is  uniformly  distributed  in  (-r,r),  we  calculate 

the  expected  value  of  y„  by  introducing  a   =  x  /2T ,    and 

setting 

1/2 

y7"  =    y   da,  ...  da    =   T  N      sgn  (xS  -  2ar)  da 

-1/2 

r 

2 


■""   '  sgn  (xS  -  x)  dx   =   i-  N[(xS  +  D  -  (F  -  xS)  ] 


-r 


=  NxS  , 


exactly  the  distance  traversed  by  the  shock,  during  the  time 


interval  t^,  =  Nx 

N 


11 


In  this  paper  we  analyze  a  modification  of  Glimm's  scheme 
where  the  exact  solution   of  the  Riemann  problem  is   eplaced 
by  an  approximate  one.   We  also  analyze  the  effect  of  replacing 
the  staggered  grid  by  a  grid  moving  with  variable,  state 
dependent,  speed. 

We  denote  by  p(x,v,w)  an  approximate  solution  at  t  =  x 
of  the  Riemann  initial  value  problem 


V   for   x  <  0 
(2.3)  u(x,0)  =  ' 

w   for   x  >  0 


for  the  equation  (1.1).   We  require  that  this  approximate 
Riemann  solver   furnish  correctly  the  integral   of  the  exact 
solution.   The  integral  of  the  exact  solution  can  be  determined 
by  integrating  the  conservation  law  (1.1)  over  the  rectangle 
(-r,r)  X  (0,t)   in  the  x,t  plane.   Denote  by  |a|  the  largest 
signal  speed  for  the  Riemann  problem  (2.3),  i.e.  assume  that 
the  exact  solution  u(x,t)  satisfies 

V   for   X  <  - i  a  1 1 


u(x,t)  = 

w   for   X  >   I  a  1 1 

Suppose  we  choose  r  >  |a|T;  then  the  result  of  the 
integration  is 

r 

(2.4a)      u(x,t)  dx  =  r(v  +  w)  +  x  f(v)  -  x  f (w)  . 

-r 

Accordingly  we  require  that 
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r 

{2.4b)      p(x,v,w)  dx   =   r(v  +  w)  +  t  f(v)  -   t  f (w) 

-r 

for  all  r  satisfying 

(2.5)  r  >  |c|t  , 

where  |c|  is  the  largest  propagation  speed  in  the  apprxximate 
Piemann  solver.   We  shall  refer  to  (2.4b)  as  oonsistenoy    of 
the    Riemann    solver  with    the    conservation    law. 

We  can  derive  an  analogous  inequality  for  the  exact 
solution  of  the  Riemann  problem  by  integrating  the  entropy 
inequality  (1.8a)  over  the  rectangle  used  before.  We  obtain 

r 

(2.6a)       U(u(x,t))  dT  <  r(u(v)  +  U(w))  +  tF(v)  -  tF(w) 

-r 

Accordingly  we  require  that  the  approximate  Riemann  solver 
satisfy 

r 

(2.6b)       U(p(x,v,w))  dx  <  r(u(v)  +  U(w))  +  tF(v)  -  tF(w) 

-r 

We  shall  call  this  property  consistency  with    the    entropy 
condition. 

Next  we  show  how  to  use  approximate  Riemann  solvers 
to  build  approximate  solutions  on  a  moving  grid.   We  denote 
by  I   .  the  j    interval  in  the  decomposition  of  the  real 
line  at  the  n    time  step.   The  decomposition  at  the  next 
time  step  is  obtained  by  moving  the  endpoints  of  I  .  with 

13 


speed  s  .  , ,  /T 
'^  nj+1/2 

To  obtain  formulas  not  overloaded  with  subscripts,  take 

the  midpoint  of  I  .as  the  origin  of  the  x  coordinate,  so  that 

■  ■  I 

I  •  =  (-r,r) 

nj     ^        r     ' 

Then  we  set 

(2-7)     I^^i^j  -    (B_,B^)  ,    B^  =  +  r  +  TSj^,/,  . 


^n+1 

/  j 

1 

1 

1 

1 

1 

1 
B_ 

1 

1 

^n: 

1 

1 

1 

-r 


We  denote  by   u^  ,  u,,  and  u„   the  constant  values  of  the 
L     M       R 

approximate  solution  on  I   .  ,  ,  I   .   and  i   .  , .   We 

n,j-l  n,D       n,j+l 

shall   build   an   approximation  to  the  solution  on 

I  ,-l  •   of  this  piecewise  constant  initial  value  problem 

out  of  the  approximate  solutions  of  the  two  Riemann  problems 

centered  at  -r,  t   and  r,  t  ;   we  assume  that  1    , -,     ■    belongs 
n         n  n+1 , ]       ^ 

to  the  domain  of  determinacy  of  I   .,ui   .ui   ...i.e, 
that  no  wave  other  than   these  two  impinge  on  I  . ,  . .  We 
define  the  approximate  solution  to  be 


,  p(x+r,u^ ,u^)  ,   B   <  x  <  0 
(2.8)     p{x,u   u  ,u  )  =  L  ^ 

'  p(x-r,Uj^,u^)  ,   0   <  X  <  B_^ 
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Integrating  the  exact  solution  of  the  conservation  law 
over  the  rectangle  In+l,j  ""    ^^n'^n+1^  gives 


J   ^(^'^n+l^  -  J   ^^^'^n^  ^^  -^  "  ^j+1/2  -  "  ^j-1/2 


=  0  , 


where   f-+i/2  ^^^  fluxes  at  the  end  points  of  I^+i^j-  If 
we  assume  that   s._^  ^,2    >  0/  it  follows  from  (2.7)  that 


^ 


j   u(x,t^)  dx  =  {2r  -  TS._^/2)u^  +  TS.^^/2^j^  . 
B 
We  call  an  approximation  p  consistent  with  the  conservation 

law  if 

(2.9)  j   p(x,Uj^,u^,u^)  dx  =  (2r-TS._^/2)Uj^+  ts.^^/2^R 
^n+l,j  +  '^fj+i/2  "  ^^j+1/2  ' 

where   f.  -.  ,-.    are  numerical    fluxes    consistent  with  the 
physical  flux,  i.e., 

(2.10)  f  .^^/2  =  f^^M'V  '    ^j-1/2  =  ^(^L'V  ' 
where  f(u,v)  satisfies 

(2.11)  f(u,u)  =  f(u)  . 

Similarly,  we  say  that  an  approximation  p  is  consistent 
with  the  entropy  inequality  if 

(2.12)  U(p(x,Uj^,Uj^,u^))  dx  <  (2r-TS^_-L/2)UM 

^-  -^  ^^j+l/2V  ^^j-1/2-  ^^j+1/2' 
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where   U^  =  U(Uj^),   U^^  =  U(Uj^)  and 

are  numerical  entropy  fluxes  consistent  with  the  physical 
entropy  flux: 

(2.14)  F(u,u)  =  F{u)  . 

Theorem  2.1.   Suppose  an  approximate  Riemann  solver 
satisfies  the  consistency  relations  (2.4b)  and  (2.6b).  Then 
the  approximate  solution  (2.8)  satisfies  the  consistency 
relations  (2.9)  and  (2.12),  with  appropriately  chosen 
numerical  fluxes. 

Proof;   Using  the  definition  (2.8) 

(2.15)  I   p(x,u^,Uj^,Uj^)  dx 

B 
-    -  ^^j+1/2 


p(y,Uj^,Uj^)  dy  +   I      p(y,Uj^,Uj^)  dy  . 

Ts .  ,  ,T  -r 

3-1/2 

We  assume  for  simplicity  that  the  speeds   s .  , /^  depend 
only  on  the  two  adjacent  states: 

(2.16)   Sj_^/2  =  s(Ul,u^)  ,    s.^^/2  =  ^(^M'^r)  ' 

where  s(u,v)  is  an  arbitrary  function.   We  now  define  two 
functions   f   and  f  : 

r 

(2.17a)  Tf-(u^,Uj^)  =  (-  r  +  TS^j^)Uj^  +  Tf^+   I   p(y,Uj^,Uj^)  dy 
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(2.17b)    Tf+(u^,Uj^)  =  ruj^  +  TSj^j^u^  +  -^^M  -   J   P(y'^M'^R^  ^^' 

-r 

Here  we  have  employed  the  obvious  abbreviations  s^..  =  s(u^,u,J, 
-^  LM       L   M 

From  the  definitions  (2.17),  and  (2.15),  we  deduce  that 


(2.18)   i   p(x,u^,Uj^,Uj^)  dx  =  (2r-TS^^)u^ 


B 


+  "^^MR^R  -^  -^^'^^L'^m)  -  ^^'"(^M'^r)- 


Lemma  2.2.   (i)   f  (u^,Uj^)  =  f'^(Uj^,Uj^)  -  fj^ 

(ii)    f"(Uj^,Uj^)   =  f'^(Uj^,Uj^)   . 

Proof ;  (i)  follows  immediately  from  the  definition 
(2.17)  upon  noting  that  p(y,u  ,u  )  =  u  .  To  prove  (ii)  we 
replace  u   ,  u*^   in  (2.17a)  by  u   ,  u   and  deduce 


(2.19)  Tf  (u^,Uj^)  -  Tf"'(u^,u^)  = 

=  (-r+Ts^j^)Uj^  -  ru^  -  Ts^^uj^  +  Tf^  -  -^^M  ^  I  P^^'^M'V  '^y 


According  to  (2.4b) ,  with  v  =  u   ,  w  =  u   ,  the  right  side 
of  (2.19)  is  zero.  This  proves  that  f~  =  f"*".   But  then  (2.18) 
is  the  desired  consistency  relation,  with  f   s  f   as  the 
numerical  flux  function.   Relation  (i)  of  Lemma  2.2  shows  that 
the  numerical  flux  is  consistent  with  the  physical  flux. 

The  proof  of  the  entropy  consistency  (2.12)  goes  the  same 
way,  except  that  instead   of  F~  =  F    we  merely  conclude  that 
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(2c20)  tF  (Uj^,Uj^)  -  tf"*'(u^,Uj^)  <  0  . 

The  analogue  of  (2.18)  holds: 
B. 


(2.21)   I  U(p(x,u^,Uj^,Uj^))  dx 


B 


=  (2r  -  TS^^)Uj^  +  TS^^U^   +  ^F-(u^,Uj^)  -  TF^u^,u^) 


Combining  (2.20)  and  (2.21)  shows  that  (2.12)  holds  with 
entropy  flux  F(u,v)  equal  to  either  F  (u,v)  or  F  (u,v) . 
This  completes  the  proof  of  Theorem  2.1. 

A  final  remark,  to  be  used  in  Section  3:   Suppose  that 
the  approximate   Riemann  solver  is  Lipschitz  continuous  in 
the  sense  that 

(2.22)  |p(x,v,w)  -  v|  <  0|v  -  w| 
for  all  x.  Then  it  follows  from  (2.8)  that 

(2.23)  |p(x,u^,u^,u^) I  <  0[|u^  -  u^l  +  \u^   -   u^ | )  . 


18 


3.     Convergence  of  the  Random  Choice  Method 

Let  I   .  denote  the  intervals  of  a  moving  grid,  as 
n,D  ^  ^ 

described  in  Section  2.   We  assume  that  the  lengths  of  these 
intervals  are  comparable,  in  the  sense  that  there  is  a  common 
length  scale  A  such  that 

(3.1)  k~-'-A  <  |l   J  <  kA 

n ,  J 

for  all  n,  j,  and  for  all  grids,  k  independent  of  A. 

Let  p(x,u^ ,u„,u^)   be  an  approximation  on  I   ,  .  to  the 
L   M   H  n+x , 3 

solution  of  (1.1)  with  initial  values  u  ,  u  ,  u  on 

L    M    K 

I   .,,1  ■  ,    1       ■  .  ■,  .      Let  T  be  the  time  at  which  we  wish  to 
calculate  the  solution,  and  denote  by  N  the  number  of  time 
steps  needed  to  reach  it,  i.e.  t   =  T.   Denote  by  a,,..., a   =  a 
N   number   contained   in  the   interval  (0,1).   The  random 
choice  solutions   v(x,t,a)   are  piecewise  constant  on  the 
intervals  of  the  grid;  let's  denote  the  value  of  v(a)  on  I 


n,D 


by  V.  =  v.(a).   These  values  are  constructed  recursively 
as  follows: 

(i)    V.  is  the  average  on  I-  .  of  the  prescribed 
initial  value  v(x,0). 

(ii)   Given  v»  for  all  £,  we  set 

,  -,  -,\  n+1     ,       n    n   n   , 

(3.2)  V.    =  p(x  ,T,v.  ,,v.,v.  ,) 

j      ^   n+1   j-1   j   j+1 

where   x  , ,  is  a  random  point  in  I  , ,  . .   Using  the  nota- 
n+1  ^         n+1,]       ^ 

tion  (2.7):   I  ,,  .  =  (B  ,  B,),   we  write 
n+1 , j      -    + 

(^•^^  ^n+1  =  %+lB-   -^  (1  -  %+l)^  • 
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We  denote  the  solution  so  constructed  by  v{x,t,a). 

We  now  impose  as  hypothesis  what  Glimm  was  able  to  prove: 
The  total  variation  in  x  of  all  approximate  solutions  v(x,t,a) 
is  uniformly  bounded  for  all  t,  all  a  and  all  N.  Under  this 
hypothesis   we  are  able  to  prove  the  following  variant  of 
Glimm' s  theorem: 

Theorem  3.1.   Let  v(x,t,a)  be  the  random  choice  solutions 
defined   as  in  (3.3)  with  the  aid  of  approximate  solutions  p 
and  random  choices  (a, ,...,a  )  =  a.  Assume  that  p  is 
Lipschitz   continuous  in  the   sense  of  (2.23),  that  p  is 
consistent  with  the  conservation  law  in  the  sense  of  (2.9), 
and  the  entropy  condition  in  the  sense  of  (2.12).  Assume 
that  the  v(x,t,a)  have  uniformly  bounded  total  variations 
with  respect  to  x.  Then  for  almost  all  choices  of  a  in  the 
sense  of  Lebesgue  measure  v(x,t,a)   converges  to  a  limit  u(x,t) 
that  has  the  prescribed  initial  value,  satisfies  the  conserva- 
tion law  in  the  sense  of  (1.3),  and  the  entropy  condition  (1.8b), 

Proof :   The  integral  form  (1.3)  of  the  conservation  law 
asserts  that 

(3.4)   -     (w^u  +  w^f)  dx  dt  -    w(x,0)  u(x,0)  dx  =  0 

for  all  C„  test  functions  w.   We  replace  w.  and  w  by  difference 
0  '^      t      X   -^ 

quotients,  and  integration  with  respect  to  t  by  summation. 
Let  us  take  for  simplicity  all  t  increments  equal: 
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(3.5a)  t  ^,  -  t   =   T 

n+1    n 

We  shall  assume  that  t  and  A  are  proportional,  i.e.  that 

(3.5b)  X   =   t/A 

is  independent  of  A.   We  define 

N     f 
(3.6)  r   =   -    I    T     \     [(D^w)  v(x,t^)  +  (D^w)  f(v(x,tj^))]  dx 


w(x, 0)  u (x, 0)  dx  , 


where 


D^w  =  [w(x,t^)  -  w(x,t^_^)]/T  , 


D  w  =  [w(x+A,t  )  -  w(x-A,t  )]/2A, 
X  n  '    n    ^       ' 

and  A  the  length  scale  introduced  in  (3.1).   Summing  by  parts 
and  shifting  the  variables  of  integration  yields 

N 

(3.7)  r  =  I  ^n 


where 


(3.8)   r^  =    w(x,t^)  [v(x,tj^_^j^)  -  v(x,t^)  +  xD^f  (v  (x,  t^)  )  J  dx 


We  take  for  the  function  v   appearing  in  (3.8)  the 

random  choice  solution  v(x,t,a);   this  makes  r   ,  and  their 

'  '  n 

sum  r,    a  function  of  a.   We  shall  show  that  when  A  is  small 
r  is  small,  except  on  a  set  of  small  measure. 
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We  need  the  following  lemma,  analogous  to  one  employed 
by  Glimm: 

Lemma  3.2. 

(i)     |r  (a) I  <  0(A)     for  all  a 
(3.9) 

(ii)  (|  r^^  da|  <_   O(A^) 

Proof;   We  break  up  r   as  the  siam 
n 

(3.10)  r   =  y  r 
where 

(3.11)  r^  .  =     f   w(x,t^) [v(x,t^,,)  -  v(x,t^)  +  TD^f (v(x,t  ))]dx 

^n+l,j 

We  deduce   from   the  definition   (3.2)  of   v(x,t  ,,)   and 

n+l 

property  (2.23)  of   p   that 


(3.12)  |v(x,t^^3^)-v(x,t^)  I  10(|v^^^-v^|  +  |v^_^  -  v^l)  . 

Since  f  is  Lipschitz  continuous, 

A 


(3.13)   x|D^f(v(x,t^))|  =  A  |f(vn^^)  -  f(v^_^) 


By  assumption  (3.1),  |£|  +  |m|  <  k.  Using  (3.12)  and  (3.13)  to 
estimate  the  right  side  of  (3.11)  gives 

|r   I  <  0(A)     y   |v5  -  v^       \. 
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Sximming  with   respect   to   j    gives 

kni    -  ^    l^ni'    -  °^^^^    ' 
J 

where  n  is  the  total    variation   of  v(x,t  ) .   Since  n  is  assumed 
bounded,  (3.9)^   follows. 

We  turn  now  to  (3.9) . .  ;  we  shall  prove  that 


(3.14)  I  j  r^  da^^3_|  <_   0(A^  , 

from  which  (3.9)..  follows. 

By  definition,  v(x,t  )  is  independent  of  c   , .  The 
dependence  of  v(x,t^_j_j^)  on  ot^^^^  is  given  by  (3.2),  (3.3). 
We  conclude  by  integrating  (3.11)  with  respect  to  o.^_^^   that 


(3.15)   I  r^^j  da^_^^  =   I     w(x,tj^)  M(x)  dx 

^n+l,j 


where 

(3.16)   M(x)  =  jY^ r       P(Y/v"_j^,v'J,v'|_^^)  dy 

'  n+1,  j  '    •'         -•     -^      -^ 
n+l,j 

-  v(x,tj^)  +  xD^f  (v(x,tj^)  )  . 
Integrate  (3.16)  with  respect  to  x  over  I 


n+1, j' 


(3.17)      I   M(x)  dx  =     J   P(y,Vj_3^,v",v^_^j_)  dy 


"'"n+l,j  "'"n+l,j 


v(x,t  )  dx  +  T        D  f  dx 
'  n  J     x 

■'•n+l,j  ^n+l,j 
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Using  (2.9)  gives  for  the  first  two  integrals  on  the 
right  the  value 


(3.18)  Tf  ._^/2  -  ^fj+1/2  ' 

where  f-.-i/2  ^^®  defined  by  (2.10).   The  third  integral  can 
be  written  as  follows: 

(3.19)  T  I   D^f  dx  =  -^  j   [f(x+A)  -  f(x-A)]  dx 

B  B 


2A 


J  +A  B_+A 

j   f(v(y,t^))  dy  -  ^    I   f(v(y,t^)) 


B_|_-A  B_-A 


^  ^j+1/2  "  ^^j-1/2 


where  f-.i/o  ^s  the  average  of  f(v(y,t  ))  over  an  interval 
of  length  2A  centered  at  the  common  boundary  of  I   ,  .  and 

■■■n+l,  j+1* 

It  follows  from  (3.1)  that  f-.-i/T  depends  only  on  v .  .  , 

kl  1  k: 

^j+1/2  =  ^(^j-k+l"-"Vk^ 

f  is  a  numerical  flux  function;  from  the  definition  of  f . , ,  ,» 
it  follows  that  this  numerical  flux  is  consistent  with  the 
physical  flux,  i.e 

f(v,...,v)  =  f(v)  . 

Setting  (3.18)  and  (3.19)  into  (3.17)  we  get 
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(3.20)         I   M(x)  dx  =   Th.^^/2  -   ^hj-1/2 


^n+l,j 


where 


^j+1/2    ^j+1/2  ~  ^j+1/2 


Since  h  is  the  difference  of  two  numerical  fluxes,  both 
consistent  with  f,  it  follows  that 

h(v,...,v)  =  0 
Since  h  is  Lipschitz  continuous,  it  follows  that 

(3.21,        |hj^,/,l  10[^^^l^    IV,  -u„|) 

Similar  estimates  hold  for  M{x) :   using   the  definition  (3.16) 

of  M  we  see  that  M(x)  =  0  when  v.  ,  =v.  ,.,  =  ...  =v.,,. 

]-k     ]-k+l  ]+k 

Using  (2.23)  we  deduce  that 

(3.22)       |M(x)|  <_   of    I       |v^    -  u  l] 

We  have  now  all  the  inequalities  needed  to  estimate  (3.15). 

Let  us  denote  the  center  of  I  , ,  .  by  x   .  and  w(x  ■  rt    )    by  w 

n+1,3  ^      n,D        n,:)'n   ^      n,j 

using  (3.20)  we  can  rewrite  (3.15)  as 


(^•23)   I  r^j  da^+i  =        [w(x,t^)  -  w^^.]  M(x)  dx 


^n+l,j 


+  Tw  .  [h .  , ,  -T  -  h .  T  /.J 
n]   D+1/2     D-l/'^ 


The  test  function  w  is  smooth,  and  by  (3.1)  the  length  of 

I   -,  .  is  0(A);  so  the  first  term  on  the  right  in  (3.23)  is  < 
n+x ,  J 

(3.24)  0(a2)     I       Iv^    -  v^l  . 

|£|<k    J^     J 
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Now  sum  (3.23)  over  j ;  since  the  sum  of  (3.24)  over  j  is  <_  0(A  )  n  r 
we  obtain   that 


(3.25)   I  I  r   da  ^, I  =  I  y  I  r  .  da  ^T  I 
'  J   n   n+1 '    '  'r  J   n]    n+1 ' 

<  0(A2n)  +  t\    I   w   [h    -  h  ,^]| 
—  ^   nj   j+-g     J  -g 

We  sum  the  second  term  by   parts,  obtaining 


(3.26)  y  (w   .  -  w   •_,!  )h.^, 
^        n,3  n,D+l   D+^ 

Using  Iw   .  -  w   -.J  <  0(A)  and  (3.21)  we  deduce  that  (3.26) 
^  '  n,]     n,3+l'  - 

2 
is  0(A)  n.   This  shows  that  (3.25)  is  bounded  by  0(A   +  AT)ri. 

Since  we  have  assumed,  see  (3.5b),   that   x  =  XA,   and  since 

n  was  assumed  bounded,  (3.9) . .  follows. 

From  Lemma  3.2   we  shall  deduce,  following  Glimm,  that 

r  =  ^  r   is  small  except  on  a  small  set.   Consider 

(3.27)  j  r^  da  =  j  [  I  r^j   da 

=  yr   da+   y    r   r   da. 
^  I   n        ii   I   n  m 
■'  nj^m  ■' 

Using  inequality  (3.9) .  we  get  that 


(3.28a)  )  I  r  ^  da  <  N  O(A^)  . 


V  f    ^ 


In  the  second  sum  suppose  n  >  m;  then  r   is  independent  of  a   -,  , 
so  that 
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r   r   da  , ,  =  r     r   da  ,  , 
n   m    n+1     m     n    n+1 


Usinq  (3.9) .  to  estimate  r   and  (3.9) . .  for  the  integral  gives 


0(A)0(A^)  ,  so 


(3.28b)         I  ^   I  r^  r^  da |  <  N^O(A^)  . 
n^m  •' 

Using  (3.28a,b)   to  estimate  (3.27)  gives 


r^  da  <  0(A) [NA  +  (NA) ^] 


Now  by  (3.5b),  NA  =  NtX  =  TA  =  0(1);  it  follows  that 

(3.29)  I  r^  da  £  0(A)  . 

Denote  by  m(e)  the  a-measure  of  the  set  where  |r(a) [  >  e . 

2  1/3 

It  follows  from  (3.29)  that  e  m(e)  £  0(A).  Taking  e  =  A  ^ 

we  deduce 

1/3 
Lemma  3.3.    |r(a)  |  f.  A  ^   for  all    except  on  a  set 

1/3 
of  measure  <_  0(A    )  . 

We  deal  with  entropy  in  an  analogous  manner.  The  integral 
form  of  the  entropy  condition   (1.8b)  is  that 


(3.30)   -  I  I  Lw^U  +  w^F]  dx  dt  - 


w(x,0)  U(x,0)  dx  <  0 


for  all  smooth,  nonnegative  test  functions  w(x,t).  The  discrete 
version  of  (3.30),  after  summing  by  parts,  is 
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N 
(3.31)  R   =    I    R 


where 


(3.32)   ^n  =  ]  w(x,t^)  [U(v(x,tj^_^^)  -  U(v(x,t^)  +  tD^F  (v  (x,  t^)  ]  dx  . 


The  analogue  of  Lemma  3.2  is 

Lemma  3.4. 

(i)    Ir  (a) I  <  0(A)   for  all   a 


(3.33) 


(ii)   I  R   da  <  O(A^) 
n     — 


The  proof  of  (i)  is  the  same  as  before;  the  proof  of  (ii) 
is  analogous;  instead  of  estimates  for  the  absolute  value  of  M 
we  have  only  upper  bounds  for  M,  and  correspondingly  only  an 
upper  bound  in  (3.33)  ..  . 

To  deduce  from  Lemma  3.4  that  R  is  bounded  from  above  by 
a  small  quantity  except  on  a  set  of  small  measure,  we  have 
to  modify  Glimm's  argument.   We  consider 

(3.34)   I  e^^  da  =  |  e^^^"  da   =  [[  Yl   e   "  daj^...da^  ; 

s  is  a  positive  paramteer  to  be  specified  later.   We  now  use 
the  fact  that  for  n  <  N-1,  R  does  not  depend  on  a„  ;  therefore 
in  the  integration  with  respect  to  a      the  factor    ~\     f  exp  sR 
is  a  constant,  and  the  integrationinvolves  only  R>j_-i- 
The  following  inequality  will  be  used: 

e'^  <  1  +  6  +  6^   for   |6  |  <  1  . 
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By  (3.33) .  ,  |r  I  £  0(A);  therefore  if  s  is  chosen  so  that 


sA  ^  0  as  A  -^  0   we  can  set  6  =  sR^_, 


and  obtain 


^^-1  2  2 


Integrating  with  respect  to  a   and  using  (3.33)  we  get 


^  ^  "-     '  1  +  sO(A^)  +  s^O(A^)  . 


We  shall  choose  s  >  1;  then 


^  e'^^-l  da   <  1  .  0(s2a2)  <  e^^^'^^) 
N  —  — 


Next  we  estimate  the  contribution  of  the  ct   ,  integration, 
and  so  on  all  the  way  to  a, .   The  result  is  the  inequality 

,-,  o^,        f   sR  ,   ^   NO(s^A^)     O(s^A) 
(3. 36)  e    da_<e         =e       , 

since  AN  is  bounded.   Let  us  denote  by  m(E)  the  measure  of  the 
set  where   R(a)  >  e;   it  follows  from  (3.36)  that 

,  ,   se  ^   cs'^A 
m(  e)  e   £  e 

The  optimal  choice  for  s  is  s  =  e/ecA,  giving 

r    \     ^       -£  /4cA 
m(e)  <_  e   ' 

Setting   e  =  A    yields 

1/3 
Lemma  3.5.   R(a)  ±  ^  for  all  a  with  the  exception  ol 

-1/3 
a  set  of  measure  <  exp(-0(A     )). 
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We  need  one  more  lemma  from  Glimm: 
Lemma  3.6.   For  all  a  and  all  m,  n. 


(3.37)      |v(t  ,a)  -  v(t  ,a)  L   <  0(|t   -  t  1)  , 
'n         m'L,  —    'm     n' 

where   I  I    denotes  the  L,  norm  with  respect  to  : 
^1  ^ 


Proof:   Integrate  (3.12)  over  I  ,,  •: 

(3.38)    I  |v(x,t^_j^j^)  -  v(x,tj^)  I  dx 

n+1  /i  ^/^/A\rin      ni.in      ni 

-'  <  0(A)   V .  , ,  -  V .   +   V   ,  -  V . 

Summing  (3.38)  over  all   j  we  obtain 


(3.39)         l^^^n+l'"^  "^^^n'^^'L   l^^^^'^ 

where  n  is  the  total  variation  of  v  at  time  t   ,  a  quantity 
that  is  assumed  to  be  uniformly   bounded.  So  (3.39)  implies 


v(t^^^,a)  -  v(t^,a)|^^  <  0(A) 


By  the  triangle  inequality  we  deduce 


[v(tj^,a)  -v(t^,a)|j^   <_0(A(n-m))  . 
Since   t  and   A  are  assumed  to  be  proportional,  we  conclude  (3.37) 
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We  now  take  a  denumerable  set  of  smooth  test  functions 
W(,  ,  1   =   1,2,...,    which  span  the  space  of  all  continuous 
functions;  we  choose  each  w„  to  be  >^  0 .   The  approximate 
solutions  v(x,t,a)  are  defined  at  discrete  times  t  .   We 
extend  their  definition  to  all  t  by  setting 

v(x,t,a)  =  v(x,t^,a)   for   t^  £  t  <  t^^^^ 

With  this  extension  the  formula  (3.6)  for  r  can  be  written  as 

(3.40)  r(a,A,£)  =  "  jj  ^^^^^z^    v(x,t,a) 

+  (D^w^)f (v(x,t,a) ]  dx  dt  -   w^(x,0)  u(x,0)  dx. 

Similarly,  we  can  write  (3.31)  as 

(3.41)  R(a,A,£)  =  -  I   [(D^W;,)  U(v(x,t,a)+  (D^w^)  F(v(x,t,a)J  dx  dt 

-  I  w^(x,0)  U(x,0)  dx  . 
It  follows  from  Lemmas  3.3  and  3.5  that 

(3.42)  |r(a,A,£)|  <  A^^^    ,  R(a,A,£)  <  b}^^ 

except  on  a  set  of  a  whose  measure  is  <  ep(A).   It  follows 

that  we  can  choose  £  as  a  function  on  A  ,  I  =   il(A),  so  that  (3.42) 

holds  for  all  £  <  £(A)  except  on  a  set  E(A)  whose  measure 

tends  to  0  as  A  tends  to  0. 
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Let  a(A)  be  a  sequence  of  points  not  in  E(A).   Since  the 
functions  v(x,t,a(A))   have  uniformly  bounded  total  variation 
in  X   and  are  uniformly  bounded,  by  compactness  we  deduce 
that  for  any  fixed  t    a  subsequence  converges   in  the  L^^ 
norm.   By  a  diagonal  argument  we  can  construct  a  subsequence 
that  is  L,  convergent  for,  say,  all  rational  values  of  t. 
Using  Lemma  3.6  we  conclude  that  this  subsequence  converges 
in  the  norm 

max  |v(t) I 
t         ^1 

to  a  limit  that  we  denote  by  u(x,t) . 

Since    a(A)  does  not  belong  to  E(A),  inequality  (3.42) 
holds  for  I   <    Z{^).      Keeping  £  fixed  and  letting  A  ->  0  in  (3.40) 
and  (3.41)  we  conclude  that  the  limit  u  satisfies  (3.4)  and 
(3.30)  for  all  w  =  w„.   Since  the  w.  span  the  space  of  all 
continuous  functions  (3.4)  follows  for  all  w;  similarly,  if 
we  choose  the  w„  so  that  every  w  >_  0   can  be  approximated  by 
linear  combinations  of  w.  with   nonnegative  coefficients, 
(3.30)  follows  for  all  w  ^  0. 

Finally  we  claim  that  if   a  e  e(A),  then  not  only  a 
subsequence  of  v(x,t,a(A))  but  the  whole  sequence  tends  to 
u(x,t) .   For  if  not,  then  by  compactness  we  could  select  a 
subsequence  that  tends  to  a   limit  u'  7^  u.   Now  u'  and  u  both 
satisfy  the  conservation    law  in  the  integral  sense  (1.3), 
and  the  entropy  condition  in  the  sense  (1.8b);  firthermore 
they  both  have  the  prescribed  initial  values  u(x,0).  Therefore 
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by  the  uniqueness  theorem  stated,  and  assumed,  at  the  end 
of  Section  1,  we  conclude  that  u'  =  u.   This  proves  that  the 
sequence  v(x,t,a(A))  converges. 

Note  that,  because  of  the  indirect  nature  of  the  conver- 
gence proof,  we  have  no  error  estimates. 

The  estimate   given  in  Lemma  3.5  for  the  size  of  the 
exceptional  a-set  is  fairly  sharp  and  shows  that  even  for 
moderately  small  values  of   A  the  probability  of  R  being  large 
is  negligible. 

We  conclude  with  a  couple  of  observations  about  deter- 
ministic schemes  over  moving  grids.   These  are  of  the  form: 

(3.43)  2rv!J"^^  =  (2r  -  TS   ,^)v!'  +  TS.   v^    +  xf .  ,^  -  xf .    , 

where   f . ^,  are  numerical   fluxes   consistent  with  physical 
flux,   see  (2.10),  (2.11).   We  further   assume  that  (3.43) 
satisfies  an  entropy  condition  of  form  (2.12): 

(3.44)  2ru(v!J"^^)  <  (2r-xs.  ,)  Viv^)    +   xs  .  , ,  U  (v!? .  ,  )  +  xF.  ,  -  xF.., 

J      —  3 — £        J         J+-1     J+-L         J   5        J+-1 

where  F .  ^j^  are  numerical  entropy  fluxes  consistent  with 
physical  entropy  flux,  see  (2.13),  (2.14).   The  estimates 
contained  in  the  proof  above  show  that  for  schemes  of  this 
type   the  analogue  of  the  convergence  Theorem  3.1  holds,  i.e. 
these  schemes  converge,  provided  the  total  variation  of  the 
approximate  solutions  are  uniformly  bounded. 
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Given  an  approximate  Riemann  solver,  we  can  construct 
a  Godunov    type    deterministic    approximation      by  setting 


/ ->  /,r-\  n+1     if     ,    n    n   n   .  , 

(3.45)     V.    ""  2f    I    P^^'^j-i'^j'^j+l^  ^^ 

^n+l,j 
where   p  is  defined  by   (2.8).   We  have  shown  in  Section  2  that 
if  the  Riemann  solver  is  consistent  with  the  conservation  laws, 
then  (2.9)  holds;  combining  this  with  (3.45)  we  obtain  (3.43). 
Suppose  further  that  the  Riemann  solver  is  consistent  with  the 
entropy  inequality;  then  according  to  Theorem  2.1  inequality 
(2.12)  holds.   Since  U  is  convex,  Jensen's  inequality  holds: 
using  (3.4  5)  we  have 

U(v^"^^)  =  u|^    J    p(x)  dxj  <  27    j    U(p(x))  dx  . 
"'"n+l,j  "'"n+l,j 

Combining  this  with  (2.12)  we  deduce  (3.44).   This  proves: 

Any  Godunov  type  scheme  on  a  moving  grid  which  is  based 
on  an  approximate  Riemann  solver   that  is  consistent  with  the 
conservation  laws  and  the  entropy  condition  converges,  provided 
that  the  total  variation  of  the  approximate  solutions  is 
uniformly  bounded. 
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4.     A  Finite  Difference  Approximation   to  the  Riemann  Problem 

In  this  section  we  shall  construct  an  approximate  Riemann 
solver  that 

(a)  is  consistent  with  the  conservation  laws; 

(b)  is  consistent  with  the  entropy  condition; 

(c)  resolves  shocks  with   absolute  sharpness. 

We  introduce  a  frame  of  reference  moving  with  velocity 
V:  y  =  X  -  Vt.   In  this  frame   the  conservation  law  (1.1) 
becomes 


(4.1)       u^  +  g(u)y  =  0  ,    g(u)  =  f  (u)  -  Vu  . 


We  approximate  the  solution  of  the  Riemann  problem  for  (4.1) 


(4.2)  u(y,0)  = 


Uj^  ,    y  >  0 


by  the  single  application  of  a  finite  difference  scheme  in 
conservation  form.  That  is,  let  g(v,w)  denote  a  numerical 
flux  consistent  with  g(u)  in  (4.1): 

(4.3)  g(u,u)  =  g(u)  . 

We  employ  a  spatial   grid  of  length  6,  smaller  than  r  but 
comparable  to  it.   The  scheme,  applied  to  initial  data   VQ(y), 
is 

(4.4a)   v(y,T)  =  VQ(y)  -  y [g (Vq (y) , Vq (y+6) )  -  g (Vq (y-6) , Vq (y) ) ] , 
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where 
(4.4b) 


y  =  T/6, 


When  v„  is  taken  of  form  (4.2)  we  obtain 


(4.5) 


^L  ' 


y  <  -6 


v(y,T)  =  ' 


^T     '      -6  <  y  <   0 


'r  ' 


'r   ' 


0  <  y  <   6 
6  <  y 


where 


(4.6a)   u^  =  u^  -  vi(g^j^  -  g^)  ,     ^r  =  ^r  "  ^(gR  "  gLR^ 


Here  we  introduce  the  abbreviations 


(4.6b)   g^  =  g(u^)  ,      g^^  =  g(Uj^)  ,      q^^  =   q(n^,M^) 


Returning  to  the  original  frame  of  reference  we  define  our 
approximation   p(x;Ut,u^)   of  the  solution  of  the  Riemann 
problem  at  time  x  by 


(4.6c)   p(x;u^,Uj^)  =  v(x-Vt,t)  =  ■ 


u  ,  x-Vt  <  -5 

J-j 

u*  ,  -5  <  x-Vt  <   0 

u*  ,  0  <  X-VT  <   6 

H 

u„  ,  6  <  x-Vt 


Lemma  4.1.   Given  r,  let  5  and  t  satisfy 


(4.7) 


6  +  Vt  <  r  ; 


then  the   approximation  (4.6c)  is  consistent  with  the 
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conservation   law   in   the    sense   of    (2.4b),    i.e. 


(4.8)  p(x;u^,u„)    dx   =    r(u^    +   u„)     -    Ttf(u„)    -    f(uj]     , 


'L'^R 


L  R' 


-r 


no  matter  what  the  numerical  flux  is,  as  long  as  (4.3)  is 
satisfied,  and  no  matter  what  V  is. 

Proof: 


-r 


p(x;u^,Uj^)  dx  =     v(x-Vt,t)  dx 

-r 


(4.9a) 


r-VT 


v(y,T)  dy  =   I  v(y,T)  dy  -  tV(Uj^-  u^^)  ; 


-r-VT  -r 

(4.7)  is  used  to  obtain  the  last  equality.  Using  (4.4a)  we  get 

r 


r  r 

(4.9b)      v(y;T)  dy  = 

-r  -r 


VQ(y)  dy  -  y   J  g(vQ(y) ,VQ(y+6) )  dy 


-r 


r-5 


r  r 

r 


=   r(u^+Uj^)  -  u 


+   U   I  g(vQ(y) ,VQ(y+6) )dy 

r 

J  g(vQ(y) ,VQ(y+6)) 


-r-6 
r 


Lr-6   -r-6 


dy 


=  r(n  _  +  u„)  -  ylg„  -  g,  ]  . 


L     R 


R    ^L^ 


Since  by  (4.1)   f  =  g  +  Vu,  setting  (4.9b)  into  (4.9a)  gives  (4.8) 
We  shall  choose  V  to  be  of  the  following  form: 


(4.10) 


V  =  Z{f^-f^)/Z{u^-   u^)  , 
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where  £  is  some  linear  functional. 

Lemma  4.2.   If  u  and  u   satisfy  the  RH  condition  (1.4)  : 

L         -ri 

f„-  f   =  S(u  -  u  )  ,  then  for  V  defined  by  (4.10)   and  g  by  (4.1), 


(4.11)         V^^l'^'r^  ^   ^    '         ^R  =  ^L  • 


Proof  is  by  direct  substitution,  and  use  of  linearity 
of  £. 

Another  consequence  of  the  definition  of  V,  of  g  and 
the  linearity  of  £  is  that  for  arbitrary  states  u,.  ,  u_, 

(4.12)  £(g^)  =  lig^)     . 

We  turn   now  to  specifying  £.   We  denote  by  U   the 
gradient  of  the  entropy  function  U.  By  definition  of  gradient 

(4.13a)  U  (u)w  =  ^  U(u  +  ew)  I 

u        d£   '         I^^Q 

We  set 


(4.13b)  i{w)    =    [U  (u_)  -  U  (uj]w  . 

UK        U    ij 


Lemma  4.4.   (a)  £(u   -  u^  )  >  0  if  u„  7^  u  . 

(b)  V  as  defined  by  (4.10)  is  uniformly  bounded. 

Proof:   We  introduce  the  linear  function  u(9): 


(4.14a)         u(e)  =  eu„  +  (i-e)u^. 

Clearly 
(4.14b)      §^   =  ^R  -  ^L  '    ^(°)  =  ^L  '    ^(^^  =  ^R' 
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Using  (4.13)  and  (4.14b)  we  can  write 


(4.15)  T 

=§0  u(u(e)) 


0      .  "uu^^R-^L^-^-R-^L^  "'' 


Being  an  entropy  function,  U  is  convex;  let  us  assume  that 
it  is  strictly  convex,  i.e.  that  U    >  cl,  c  >  0.  It  then 
follows  from  (4.15)  that 

(4.16)  £(u^-  u^)  >  c|u^-  u^l^. 

Lemma  4.4  follows  from  this  lower  bound. 

Our  choice  for  numerical  flux  is  of  Lax-Wendroff  type: 

(4.17a)         g^^  =   9q(u^,u^)    -    B(Uj^-u^) 

where  g,.  is  van  Leer's  modification  of  Richtmyer's  two-step 

scheme : 

u  +u„ 
(4.1 /b)  w  =  -^2-^  -  ^  [9^-9^]     , 

(4.17c)     gQ(Uj^,Uj^)  =  g(w)  -  g( — ^ )  +  — 2 ' 

see  (4.6b)  for  the  abbreviations. 

The  following  estimate  is  an  easy  consequence  of  the 
definition   g„  ,  the  Lipschitz  continuity  of  f  and  the  uniform 
boundedness  of  V: 

(4.18)        Igo(uL.u^)  -  gj  10(|uj^-  uj). 


39 


The  coefficient  B  in  (4.17a)  is  an  artificial  viscosity 
coefficient;  we  take  it  to  be  scalar.  Of  course  it  must  satisfy 

(4.19)  0  <  B(Uj^,Uj^) 

Equations  (4.6),  (4.10),  (4.13)  and  (4.17)  define  our 
approximate  Riemann  solver.   The  rest  of  this  section  is 
devoted  to  specifying  6(u,v)  so  that 

(i)   Shocks  are  resolved  sharply. 

(ii)  The  entropy  condition  (2.6b)  is  satisfied 

Lemma  4.5.   The  approximate  Riemann  solver  resolves 
shocks  sharply  iff 

(4.20)  6(u^,u^)  =  0 
wheenver  u   and  u   are  connected  by  a  shock. 

Ij        K 

Proof:   If  u  ,  u    are  connected  by  a  shock  traveling 
with  speed  S,  then  the  RH  condition  (1.4)  is  satisfied. 
According  to  Lemma  4.2,  V  =  S  and  g„  =  g^.  .   Setting  this  in 

K         Li 

(4.17b)  gives   w  =  (Ut+Ut,)/2  ;  setting  this  in  (4.17c)  gives 

(4.21a)  ^O^^L'V  =  ^L  =  ^R 

Setting  this  in  (4.17a)  gives 

^LR  =  ^L  -  ^("r-  ^L^  =  5r  -   ^^^R-  ^L^  • 
Substituting  this  into  (4.6a)  gives 


(4.21b) 


u^  =  u^  +  yB(Uj^  -  u^) 


^R  =  ^R  -  ^^(^R  -  ^L^ 
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Sharp  resolution  means  that  the  random  choice  method  intro- 
duces no  intermediate  values.  This  is  the  case  only  if 

*         * 
u   =  u   ,  u   =  u  .   This  proves  Lemma  4.5. 

J_i       Li       K      R 

We  now  turn  to  the  entropy  condition  (2.6b) : 


(4.22)   J  U(p(x,Uj^,Uj^)  )  dx  -  r[U(Uj^)  +  U(Uj^)]  +  tF  (u^^) -tF  (u^^)  <  0, 

-r 

Using  the  definition  (4.6c)  of  p  we  obtain 


U(p(x,Uj^,u^)  )  dx 


-r 


(4.23a) 


=  (r-6+VT)U(uJ+  6U(u*)  +  6U(u*)  +  (r-6-VT)  U(u^) 

ij  J-j  K  K 


=  r[U(u^)  +  U(Uj^)]  +  6[U(Uj^)  -  U(L)  +  U(Uj^)  -  U(Uj^)] 


+  TV[U(u^)  -  U(u^)  ]  . 
Let  us  denote  the  entropy  flux  in  the  moving  frame  by  G: 

(4.23b)  G(u)  =  F(u)  -  VU(u)  . 

Setting  (4.23)  into  (4.22)  we  obtain: 


(4.24)   6  [U(Uj^)-U(Uj^)+  U(Uj^)-U(Uj^)]  +  t[G(Uj^)  -  G(Uj^)]  £  0 


Our  aim  is  to  choose  3  so  that  (4.24)  holds,  while 
retaining  the  restrictions  previously  placed  on  3.   Denoting 
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the  subscripts  L  or  R  by  the  common  symbol  P,  we  write  the 
linear  approximation  to  U: 


(4.25)  U(u*)  =  U(Up)  +  U^(Up) (Up  -  Up)  +  Q  . 
According  to  (4.6a) 

u*  -  u^  =  -y(gLR-  gL^  '  ^r  "  "r  =  ^(^lr"  ^r^ 

We  set  this  into  (4.25);  using  the  abbreviation 

(4.26)  U^(Up)  =  Vp 

we  can   rewrite  (4.25)  as 

(4.27)  U(u*)  -  U(Up)  =  +  -Vp(g^j^-  g^)  +  i^^Q^ 
where 

(4.2a)  iQpl  1  o(|g^j^-  gpl^) 

Setting  (4.27)  into  (4.24)  gives,  after  division  by  t  (recall 
that  y  =  6t  )  : 

(4.29)  y(Q^+  Q^)  -  v^(g^j^-  g^^)  +  v^(g^^-  g^^)  +  G^-  G^  <  0  . 
Next  we  set  into  (4.29)  the  definition  (4.17a)  of  9x0- 

(4.30)  <,^^  =   g^  -  B(u^-  u^)  . 

After  a  little  rearrangement  of  some  of  the  terms  we 
get  the  following  inequality: 


(4.31)  yH-j^  +  H2  +  H3  <  3H^ 
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where 


(4.32)  ^  ^^ 


«4  =  ^^R-  ^L^  (^R-  ^L' 


We  estimate  the  first  three  quantities   from  above,  the 
last  one  from  below: 

Lemma  4.6.   Denote  by  K  the  quantity 


{l^o  -  ^l'  -"   1^0  -  ^r|}i^r  -  ^I 


(4.33)     K  -  ng^  -  gj  +  Ig^  "  9 


Then,  with  suitable  constants  C.  , 


H^  <  C^K  +  C^B^IUj^  -  u^\^ 


"2  1  C^K 

(4.34) 

«3  iS'^R  -  ^lI' 


H4  1  c\u^   -    u^|2 


Proof;   Setting  the  definition  (4.30)  of  g    into 
inequality  (4.28)  shows  that 


"1  i°1i9o-  9lI^  *  lgo-  %l^  *  ^^ 


l"L-  -rI'} 
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Using  inequality  (4.1b)  gives  (4.34)  j^.   Inequality  (4.34)2 
follows  from  the  Schwarz   inequality  applied  to  H^   as 
defined  in  (4.32)2  '  using  [  v^^-  v^^  |  <  0  (  |  u^^-  u^  |  )  . 
To  estimate  H   we  start  with  the  definition  (4.23b)  of  G: 


S  -  ^L  =  ^R  -  ^L  -  ^("r  -  "l'  • 


Employing  the  linear  function  u(e)  defined  in  (4.14a)  we 
can  write  this  as 


t^u  -  ^"u^  ("r  -  "l^  ^'    • 


Using  relation  (1.6)  between  F^  and  U^  ,  and  then  the 
definition  (4.1)  of  g,  this  can  be  rewritten  as 


U  If   -   V]  (Up  -  u  )  d( 
u   u         R     L 


U  g  (u„  -  u  )  d( 
u^u   R     L 


u  de 


Integration  by  parts  turns  this  into 


(4.35a)     G  -G,  =  Vp(g  -g,)  -    U  ^,  (u  -u,  )  •  (g  -g,)  d 


R   L     R'^'R  ^L 


uu   R   L 


Setting  this  into  the  definition  (4.32) ^  of  H^  gives 


"3  "    2 


'R   ^L' 


(go-  gr)  -     u  n^^p-  ^T)-(g  -  ^t.) 


uu   R    L 


Inequality  (4.34)^  is  an  immediate  consequence. 
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Finally,  setting  (4.26)  into  (4.13b)  we  conclude  that 

(4.35b)  H^  =    Ku^   -   u^)     ; 

therefore  (4.34)   is  just  inequality  (4.16). 
We  now  define  B  by 

(4.36)  H^6  =  2C  K  +  2  H^  , 

where   H-,  =  max  (H^,0)  . 

Lemma  4.7.   With  B  defined  as  in  (4.36),  the  entropy 
inequality  (4.31)  holds,  provided  that  C  is  taken  large 
enough  and   y  small  enough. 

Proof:   By  inequality  (4.18) 

K  iCju^  -  u^|2  . 

Using  this,  and  inequalities  (4.34) 2,  (4.34)^   we  conclude 
from   (4.36)  that 

6  <  (2CC^  +  2C^)/c    =   C5  . 

Setting  this    into  (4.34)   yields 

Using  this  and  (4.34) _  gives  an  upper  bound  for  the  left 
side  of  (4.31).   Thus  (4.31)  follows  from 


(4.37a)   {\iC^    +   C2)K  +  mC^C^BIu^^  -  u^^  |  ^  +  H3   £   6  H^ 
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Clearly   H   <_  H^  ;  so,  by(4.36), 

(4.37b)  {VC^+   C^)K   +  H3  <  I  BH^ 


provided  that 


{4.38a)  pC,  +    Cj    <_  C 


On  the  other  hand   using  (4.34) .  we  conclude  that 
(4.37c)  yCj^C^BlUj^  -  u^l^  <  I  3H4 

provided  that 


(4.38b)  yC^C^  -k  ^ 

Clearly,   (4.38a)  is  satisfied  if  C  is  large  enough,  and 

(4,38b)  if  y  is   small  enough.   Adding  (4.37b)  to  (4.37c) 

yields  the  desired  inequality  (4.37a),  and  completes  the 
proof  of  Lemma  4.7. 

Remark:   The  entropy  inequality  provides  an  a  priori 
estimate  for  solutions,  and  therefore,  guarantees  stability. 
Condition  (4.38b)  is  a  kind  of  CFL  condition;  this  can  be 
made  more  transparent  in  specific  cases,  see  Sections  5  and  6. 

We  verify  now  that  our  choice  of  B  in    (4.36)  is 
consistent  with  previous  conditions  imposed. 

(a)  Clearly,   6^0,  satisfying  (4.19). 

(b)  Suppose  u   and  u   are  connected  by  a  shock. 

Then  the  RH  condition  holds,  so  according  to  (4.21a)  of 

Lemma  4.5,  q^  =  g^  =  g^,.   Setting  this  into  the  definition 
^0     L     R 
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(4.33)  of  K  and  (4.32)  of  H^   we  get 

(4.39)  K  =  0  ,    E^    =   G^   -    G^    . 

According  to  Lemma  4.2,  if  the  RH  conditions  hold  then  V   -   S, 
setting  this  into  the  definition  (4.23b)  of  G  we  get 

S-  ^L  =  ^R-  ^L  -  ^(Ur  -  "l) 

According  to  the  entropy  condition  (1.8c),  this  quantity 
is  <^  0.   Combining  this  with  (4.39)  we  see  that  H-.  <_  0, 
so  H-,  =  0;  setting  this  into  (4.36)  we  get  3  =  0,  as 
required  in  (4.20). 
We  summarize: 

Theorem  4.8.   The  approximate  Riemann  solver  defined 
by  (4.6),  (4.10),  (4.13),  (4.17),  (4.33),  (4.34)  and  (4.36) 
is  consistent  with  the  conservation  laws,  the  entropy  condi- 
tion, and  resolves  shocks  sharply,  provided  that  the  CFL 
conditions  (4.7),  (4.38b)  are  satisfied. 

In  conclusion   we  note  that  our  scheme  has  the  following 
monotonicity  property: 

Theorem  4.9.   The  approximate  Riemann  solver  is 
monotonic  in  the  sense  that 

ii(Uj^)  <^  ^(^l)  1  ^^^R^  -  ^^^R* 

Proof:   Since  £,  is  linear,  these  inequalities  can  be 
rewritten  as 


(4.40)   0  <  i{u^   -    u^)  <  £(Uj^  -  u^)  <  £(Uj^  -  u^)     . 
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Using  the  definition  (4.6a)  of  u   ,  u   ,  and  definition 

J_i  X\ 

(4.17a)    of   g   we   have 


u^   -   u^   =   -    y(gQ   -   g^)    +   v)B(Uj^  -   u^) 


Setting   this    into    (4.40)    gives 

0    <    v&Ziu^   -   u    )    -    vUgf.    -   gj 


<    (l-y6)£(Uj^-   xi^)    +   M£(gQ-   gj^) 


<    £(u^   -   u^)     . 

The  first  and  third  inequality   are,  on  account  of  (4.12), 
the  same  and  can  be  written  in  the  symmetric  form 

9  ■'"'3 
(4.41a)         £(gQ  -  ^.^   ^)    <  BHu^   "  "l  ^  * 

The  middle  inequality  is 

g   g 
(4o41b)   -2y£(gQ  -   ^~  ^  )   1   (1  -  2y3)  Ku^-    u^)     . 

Using  definitions  (4.32) „  and  (4.35b)   we  can  write  these  as 

"2  -  ^"4  '    "^^"2  1  ^^  "  2yB)H^  . 

Since  by  (4.31),  |h„|  <_  BH.  ,  the  first  inequality  holds 
without  any  further  condiition,  and  the  second  under  the 
CFL  type  condition 

4yB  <  1. 
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5.     Scalar  Conservation  Laws. 

In  this  section  we  consider  scalar  conservation  laws; 

(5.1)  u^  +  f^  =  u^  -f  au^  =  0  ,    a(u)  =  ||  , 

u,  f  and  a  scalar  quantity. 

In  the  scalar  case  the  RH  condition  (1.4)  can  be 
satisfied  for  any  pair  u   and   u   .   It  then  follows  that 

f_  -  f, 

(5.2)  V  =  S  =  ^jS__^ 

R     L 

and  from  (4.21a)  that 

(5.3)  go  =  q^=   g^  . 

In  the  scalar  case  any  convex  function  U(u)  is  an 
entropy  function;  for  in  this  case  equation  (1.6)  can  be 
satisfied  for  any  U  by  an  appropriate  F.   To  analyze  the 
entropy  conditions  we  write  (4.35a),  using  (5.3)   as 
follows : 


(^•^)    S  -  ^L  =  -   "uu^^)  lg(u)-gj  <Ur-u^J  de  , 


where   u  =  uO)  =  Bu   +  (l-0)u^. 

R  •I-' 

Suppose  u   and  u   are  the  values  on  the  left  and  right 
L        R 

of  a  discontinuity  in  a  weak  solution  of  (5.1)  that  satisfies 
the  entropy  condition  (1.8c)  for  the  entropy  U.  This  condition 
can  be  written  as 


S  -  ^L  i  0  • 
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Clearly,  it  follows  from  (5.4)   that   this   is   satisfied 
for  all  convex  functions  U  iff 

g(u)  -  g 

(5.5)    ^ >  0  for  all  u  between  u_  and  u_. 

u^  -  u^   -  L       R 

This  condition  is  due  to  Oleinik,  who  derived  it  in  a 
somewhat  different  way.   It  is  well  known  that  solutions 
of  (5.1)  in  the  integral  sense  whose  discontinuities  satisfy 
(5.5)  are  uniquely  determined  by  their  initial  values. 

We  turn  now  to  our  approximate  Riemann  solver;  the 

formulas  are  greatly  simplified,  on  account  of  (5.3): 

*       * 
(i)   u   and  u   are  given  by  (4.21b) : 

Li  K 


(5.6)   ^L  =  ^L  +  yli(Uj^-  u^)  ,    Uj^  =  Uj^  -  y6(Uj^-  u^)  . 
(ii)  K,  as  defined  by  (4.33)   and  H   ,  as  defined 


by  (4.32)  T  ,   are  equal  to 


(5.7)  K  =  0  ,    H^  =  Gj^  -  G^  . 

With  these  values  the  choice  of  6  suggested  by  formula  (4.36)  is 

(5.8)  2R"^ 

where,  using  (5.4)  and  (4.15),  (4.35b), 


,  "uu^^-^L^^^R-^l)  ^' 
"3 

(5.9)  R  =  rr^  = 


^4  r  2 

^        U   (u^-  u^)   dG 
J   uu   R    L 

The  entropy  condition  (4.31)  will  be  satisfied  as  long 

as  3  >  (5.8).   Since  we  wish  to  satisfy  the  entropy  condition 
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for  all  entropies,  we  shall  choose  B  as  the  maximum  of 
2R   among  all  entropies;  using  (5.9)  we  get 

g(u)  -  g. 

(5.10)  B   =  -  2  min 


where  u  varies  between  u   and  u  . 

J_i         K 

Note  that,  since  the  quantity  to  be  minimized  is  zero 
at  u  =  u   ,  the  value  of  the  minimum  is  <_   0;  this  shows  that 
3  >^  0,  as  required.   Note  also  that  3  =  0  if,  and  only  if  the 
discontinuity  u   ,  u   is  a  shock,  i.e.  it  satisfies  Oleinik's 
entropy  condition  (5.5).   In  fact,  the  value  of  3  is  exactly 
twice  the  maximum  amount  by  which  the  Oleinik  condition 
is  violated.   In  actual  applications  we  did  not  bother  to 
evaluate  the  minimum  in  (5.10)  but  repacced  it  by  the  value 
of  (g  -  g, )/(u^-  u^ )   at,  say,  the  midpoint  u  -    (u  +  u  )/2. 

Note  further  that  it  follows  from  (5.6)  that  if  y  is 

chosen  so  that  yB  f^  1/2,  then  the  sequence  of  approximate 

*     * 
values   u   ,  u   ,  u   ,  u   is  monotonic.   We  deduce  from 

Li      J-i      K      R 

this  that  the  total  variation  of  approximate  solutions  at 

time  t  .,  is  <  the  total  variation  at  time  t  .  It  follows 
n+1     —  n 

that  the  total  variation  of  all  approximate  solutions  is 
less  than  the  total  variation  of  the  initial  values,  for 
all  time  and  for  all  choices  of  the  random  parameters  a. 
Using  Theorem  3.1   we  deduce  the  convergence  of  our  random 
choice  method  for  almost  all  choices  to  a  solution  of  the 
initial  value  problem  that  satisfies  Oleinik's  condition, 
provided  that  the  prescribed  initial  values  have  bounded  total 
variation. 
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Next  we  present  numerical  solutions  of  the  following 
initial  value  p  roblem: 

(5.11)       u^  +  (u^/2)   =  0  ,    u(x,0)  =  2  +  sin  x  . 

The  exact  solution  of  this  problem  is  smooth  for  t  <  1. 
At  time  t  =  1  a  shock  is  beginning  to  form  which  at  first 
grows  in  strength,  then  decays. 

In  Tables  I  and  II  we  list  along  side  each  other  the 
exact  solution,  the  random  choice  solution  on  a  staggered 
mesh,  and  the  approximation  obtained  by  the  Lax-Friedrichs 
difference  scheme.   In  both  approximations  X  is  determined 
by  the  CFL  condition,  i.e.,  at  each  time  level  At  is 
calculated  by 

-7 max  |a(u  )  I  =  1  ,    Ax  =  u/lO  . 

Ax         I  ^   '  I      '  / 

Table  I  shows  solutions  at  t  =  0.628  (6  full  time  steps) 
at  which  time  the  exact  solution  is  still  smooth.  The  relative 
ilp-error  of  the  random  choice  solution  is  about  9%;  the 
corresponding  error  of  the  Lax-Friedrichs  solution  is  about 
4%,  which  is  comparable  in  size.   Table  II  shows  solutions 
at  t  =  5.725  (50  full  time  steps),  at  which  time  the  shock  has 
already  started  to  decay.   We  observe  that  the  shock  in  the 
random  choice  finite  difference  solution  is  completely  sharp, 
and  the  values  on  the  two  sides  of  the  shock  are  accurately 
reproduced;  the  location  of  the  shock  differs  by  r  from  its 
exact  location.   Observe  also  that   the  rarefaction  region  in 
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the  approximate  solution  has  a  step-like  structure  containing 
less  infprmation  than  is  present  in   the  exact  solution. 

The  calculations  in  the  able  made  use  of  a  recipe  for 
choosing  the  random  parameters  a   recommended  by  Chorin  in  [  2 ] 
We  found  this  choice  very  helpful  for  improving  the  accuracy 
with  which  the  shock  is  propagated,  but  no  help  at  all  in 
eliminating  the  step-like  structure  in  the  rarefaction  wave. 
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Table    I.       Solutions   of    (5.11)    at    t=    0.62! 


j 

1 
X 

Exact 
Solution 

Random 

Choice 

Lax-Fried 

richs 

Solution 

Error 

Solution 

Error 

2 

-  2.827 

.300E+01 

.300E+01 

.142E-13 

.299E+01 

-.621E-02 

3 

-  2.513 

.294E+01 

.295E+01 

.127E-01 

.286E+01 

-.777E-01 

4 

-  2.199 

.267E+01 

.200E+01 

-.671E+00 

.246E+01 

-.212E+00 

5 

-  1.885 

.200E+01 

.169E+01 

-.309E+00 

.192E+01 

-.833E-01 

6 

-  1.571 

.133E+01 

.119E+01 

-.138E+00 

.148E+01 

.146E+00 

7 

-  1.257 

.106E+01 

.105E+01 

-.127E-01 

.122E+01 

.160E+00 

8 

-   .942 

.lOOE+01 

.lOOE+01 

.623E-03 

.112E+01 

.123E+00 

9 

-   .628 

.104E+01 

.105E+01 

.119E-01 

.113E+01 

.861E-01 

10 

-   .314 

.114E+01 

.120E+01 

.549E-01 

.120E+01 

.557E-01 

11 

-   .000 

.128E+01 

.141E+01 

.126E+00 

.131E+01 

.310E-01 

12 

.314 

.144E+01 

.143E+01 

-.121E-01 

.145E+01 

.106E-01 

13 

.628 

.162E+01 

.170E+01 

.819E-01 

.161E+01 

-.635E-02 

14 

.942 

.181E+01 

.200E+01 

.196E+00 

.179E+01 

-.202E-01 

15 

1.257 

.200E+01 

.230E+01 

.296E+00 

.197E+01 

-.309E-01 

16 

1.571 

.219E+01 

.233E+01 

.134E+00 

.215E+01 

-.381E-01 

17 

1.885 

.238E+01 

.233E+01 

-.507E-01 

.234E+01 

-.412E-01 

18 

2.199 

.256E+01 

.260E+01 

.388E-01 

.252E+01 

-.395E-01 

19 

2.513 

.272E+01 

.281E+01 

.901E-01 

.269E+01 

-.325E-01 

20 

2.327 

.286E+01 

.281E+01 

-.434E-01 

.284E+01 

-.201E-01 

21 

3.142 

.296E+01 

.295E+01 

-.700E-02 

.925E+01 

-.540E-02 
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Table  II.   Solutions  of  (5.11)  at  t  =  5.725 


j 

X 

Exact 
Solution 

Random 

Choice 

Lax-Fried 

richs 

Solution 

Error 

Solution 

Error 

2 

-  2.827 

.175E+01 

.174E+01 

-.682E-02 

.173E+01 

-.156E-01 

3 

-  2.513 

.179E+01 

.179E+01 

-.388E-02 

.177E+01 

-.192E-01 

4 

-  2.199 

.184E+01 

.179E+01 

-.503E-G1 

.182E+01 

-.219E-01 

5 

-  1,885 

.189E+01 

.181E+01 

-.743E-01 

.i86E+01 

-.242E-01 

6 

-  1.571 

.193E+01 

.199E+01 

.566E-01 

.191E+01 

-.262E-01 

7 

-  1.257 

.198E+01 

.205E+01 

.697E-01 

.195E+01 

-.280E-01 

8 

-   .942 

.203E+01 

.206E+01 

.370E-01 

.200E+01 

-.297E-01 

9 

-   .628 

.207E+01 

.208E+01 

.102E-01 

.204E+01 

-.312E-01 

10 

-   .314 

.212E+01 

.211E+01 

-.129E-01 

.209E+01 

-.327E-01 

11 

-   .000 

.217E+01 

.228E+01 

.113E+00 

.213E+01 

-.340E-01 

12 

.314 

.221E+01 

.228E+01 

.682E-01 

.218E+01 

-.351E-01 

13 

.628 

.226E+01 

.233E+01 

.681E-01 

.222E+01 

-.362E-01 

14 

.942 

.231E+01 

.236E+01 

.521E-01 

.227E+01 

-.372E-01 

15 

1.257 

.235E+01 

.239E+01 

.392E-01 

.231E+01 

-.382E-01 

16 

1.571 

.240E+01 

.241E+01 

.157E-01 

.236E+01 

-.410E-01 

17 

1.885 

.244E+01 

.150E+01 

-.945E+00 

.232E+01 

-.123E+00 

18 

2.199 

.156E+01 

.152E+01 

-.379E-01 

.197E+01 

.413E+00 

19 

2.513 

.161E+01 

.168E+01 

.761E-01 

ol7lE+01 

.106E+00 

20 

2.827 

.165E+01 

.174E+01 

.830E-01 

.167E+01 

.142E-01 

21 

3.142 
1 

.170E+01 

.174E+01 

.371E-01 

.169E+01 

-.841E-02 
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